Screening of Spherical Colloids beyond Mean Field — 
A Local Density Functional Approach 
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We study the counterion distribution around a spherical macroion and its osmotic pressure in 
the framework of the recently developed Debye-Hiickel-Hole-Cavity (DHHC) theory. This is a 
local density functional approach which incorporates correlations into Poisson-Boltzmann theory by 
adding a free energy correction based on the One Component Plasma. We compare the predictions 
for ion distribution and osmotic pressure obtained by the full theory and by its zero temperature limit 
with Monte Carlo simulations. They agree excellently for weakly developed correlations and give the 
correct trend for stronger ones. In all investigated cases the DHHC theory and its computationally 
simpler zero temperature limit yield better results than the Poisson-Boltzmann theory. 
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INTRODUCTION 



The screening of charged macromolecules in an elec- 
trolyte solution is a long standing problem which has 
prompted many attempts aiming at a theoretical expla- 
nation. In their pioneering work Gouy [| and Chapman 
used what is now referred to as Poisson-Boltzmann 
(PB) theory as the basis for a mean field treatment of 
the electrical double layer. This approach found its cul- 
mination about thirty years later in the famous DLVO 
theory of charged colloids bJQ. The major flaw of these 
mean field approaches is their neglect of correlations be- 
tween the ions. The first attempt to work out such 
correlations for homogeneous electrolytes are due to De- 
bye and Huckel whose work remarkably (and at first 
glance confusingly) is also based on (linearized) Poisson- 
Boltzmann theory. In the inhomogeneous case integral 
equation theories 0, Q, S an d recently field theories 
10] have become very popular in calculating correlation 
corrections to mean field double layers. However, in or- 
der to make progress and calculate physical quantities, 
approximations have to be made which, unfortunately, 
instead of clarifying the physics sometimes tend to ob- 
scure it. Moreover, since in some of these methods, the 
free energy is not defined in a unique way, it becomes 
impossible to determine the specific role played by each 
source of correlations in the system. 
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It would therefore be desirable to have a theoretical 
framework which retains the simplicity of the early at- 
tempts, but also accommodates correlation effects. This 
is the case for density functional theories. It is possible 
to rigorously rewrite the partition function of, say,a sys- 
tem of charged colloids, as a density functional [Tj, in 
which the contribution beyond mean field is seen to be 
expressible as an additive correlation correction to the 
free energy density, whose functional form is of course 
unknown and for which one has to make a reasonable 
ansatz. The spirit is very similar to the fundamental 
problem of integral equations, where one also has to make 
an educated guess (namely, the closure relation), but in 
the functional case the ansatz involves a free energy den- 
sity rather than a relation between two- and three-point 
functions. It thus relies on a different kind of intuition 
and thus permits complementary insight. 

One suggestion for such a functional correction has 
been made by Nordholm [l^. It relies on a Debye- 
Hiickel treatment of the One Component Plasma (OCP) 
[hH Ull Il5| |. in which the short-distance failure of lin- 
earization is cleverly overcome by postulating a correla- 
tion hole. Since beyond a certain density the resulting 
OCP free energy density is a concave function of den- 
sity, this favors the development of inhomogeneities. In 
the pure OCP these are balanced by the homogeneously 
charged background. However, if one uses the OCP free 
energy density as a correlation correction to the mean 
field functional describing the double layer at a charged 
surface, one has all the charge opposite to the counterions 
located on that surface, rather than homogeneously dis- 
tributed as a stabilizing background. The consequence is 
that the double layer becomes unstable and all ions col- 
lapse onto the surface, an effect which has been termed 
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"structuring catastrophe" 0, . 

To circumvent this instability without losing the phys- 
ical transparency of a local functional, we recently pro- 
posed the Debye-Hiickel-Hole-Cavity (DHHC) theory 
|l8| . in which we suggested a convex correlation func- 
tional. This was achieved by excluding the homogeneous 
background from a region of radius a around the central 
ion during the Debye charging process. For counterions 
with size we identified a tentatively as the ion diameter. 
We then applied our theory to the screening of a charged 
rod by its counterions. Comparisons of the ionic charge 
distribution obtained showed a very good agreement with 
the simulations for both monovalent and trivalent coun- 
terions. 

In this paper we test our theory for a different geom- 
etry: charged spherical colloids with point-like counte- 
rions. In general, colloidal systems exhibit a rich phase 
behavior. The particles can agglomerate at high densi- 
ties, generally an irreversible process, but they may also 
show a reversible liquid-vapor phase separation similar 
to the one present in simple molecular liquids. In order 
to prevent them from simply falling out of solution, one 
needs some kind of repulsion between the particles. Intro- 
ducing charged groups at the surface of the colloid is one 
way to do that. The large gain in entropy following the 
dissociation of a vast number of counterions into solution 
stabilizes the system, because an aggregation of colloids 
into a small sub-volume would - for reasons of global 
charge neutrality - also require the counterions to occupy 
this small volume and thereby give up much entropy. Of 
course, the final state of the system is always a balance 
between energy and entropy, and if electrostatic interac- 
tions are strong, they will ultimately overcome ent ropy 
and lead to aggregation of the colloids pH El I2I I2I l23f . 
The resulting phenomenon of "like charge attraction" has 
received much attention, but it is of course only mysteri- 
ous if one forgets that the entire system is neutral. Ad- 
mittedly, confusion persists about whether such a phase 
separation could also happen within mean field theory. 
Even though rigorous proofs exist that PB theory will not 
permit attraction between like charged macroions under 
reasonably general circumstances j24[ IH |2(| , and that 
in a cell model treatment the compressibility will be pos- 
itive |27| . it has been claimed that an expansion of the 
free energy of a charged colloidal suspension into zero- 
, one-, two- etc. body terms will contain configuration 
independent volume terms, which may drive a phase sep- 
aration even though the pair terms are purely repulsive 
pslEflHsnj l . Since unfortunately all these derivations rely 
on a linearization of PB theory l which might render the 
findings as artifacts |3j |33, yj] , the issue appears to 
be open yet. 

All these phenomena ultimately depend on the screen- 
ing produced by the ionic cloud, which in turn depends 
on the geometry of the system. In this regard, a charged 
spherical colloid differs from a charged rod in two funda- 
mental ways: the electrostatic potential and the spatial 
extension. The logarithmic potential present in the case 



of charged cylinders leads to the phenomenon known as 
Manning condensation j^, Elj . If the line charge den- 
sity exceeds a critical threshold, a certain fraction will 
remain loosely associated with the rod, even at infinite 
dilution, and renormalize the rod charge. A quantita- 
tive PB treatment of this provides a unique criterion for 
defining the effective charge of the system, even at finite 
densities [37113^. 

The situation is different for charged spherical colloids, 
which lose all their counterions in the limit of infinite 
dilution; thus, the colloidal charge does not get renor- 
malized. Still, on often talks about effective charges, 
which mimic the stronger condensation of nonlinear the- 
ory within a linearized treatment [HE1 El El El 
That, however, is clearly not a physical but rather a for- 
mal renormalization, necessitated by the simplified linear 
treatment, and is thus a different story. 

Another important difference between the spherical 
and the cylindrical symmetry lies in the spatial extend. 
If a charged rod is infinitely long (as is usually assumed in 
theoretical treatments), the number of counterions at any 
given distance from the rod is always infinite. In contrast, 
for a charged spherical colloid the number of counterions 
at any distance is always finite, since of course there is 
no direction along which the colloid is infinite. Hence, 
fluctuations of the radial charge density are more likely 
to be important in the spherical case. 

The systems we will consider here are strongly charged 
colloids with point-like ions of some specific valence and 
no added salt inside a spherical cell. Since all the parti- 
cles are limited to be within one cell, correlations between 
different macroions and between microions belonging to 
different cells are not present. In our treatment we will 
thus exclusively focus on questions regarding the descrip- 
tion of a single double layer. Furthermore, for point-like 
ions the interpretation of our cutoff parameter, a, can 
obviously no longer be the particle diameter. We will in- 
troduce an alternative prescription for a, based again on 
local density considerations and keeping in mind that its 
entire purpose is to prevent the functional from becoming 
unstable. 

We also derive an approximated version of our corre- 
lation functional, namely, its zero temperature limit. It 
has the huge advantage that it can be calculated ana- 
lytically, while still predicting ion profiles quite close to 
the full DHHC expression for a wide range of parame- 
ters. It also demonstrates the spirit of our stabilization 
correction very directly. 

Finally, we compare our predictions for ion profiles 
with Monte Carlo simulations, in which we independently 
vary valence v and plasma parameter T2d = \J ttctI^v 3 , 
where a is the density of surface charges and £b is 
the Bjerrum length. It has been shown that beyond 
T2d — 2.26 the force-distance curves between charged 
plates cease to be monotonic, and beyond T2d — 2.45 at- 
tractions set in |45| . These effects result from correlations 
between different double layers (like, for instance, ion in- 
terlocking El E3)j which we cannot account for, and 
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it has in fact been shown that they cannot be described 
within a local density functional theory with a convex 
correlation correction j2||. However, for the description 
of a single double layer the regime of applicability of our 
theory is larger, even though it clearly must fail for too 
high coupling. 

The paper is organized as follows. In Sec. II the DHHC 
correlation functional is revisited and its zero tempera- 
ture limit is introduced. It is then applied as a local cor- 
relation correction to the problem of screening of charged 
colloids in Sec. III. The case of point-like ions is discussed 
in detail and the new expression for a is proposed. Tech- 
nical details of the simulations are described in Sec. IV. 
The results of the simulations, full theory and zero tem- 
perature limit are compared in Sec. V, and we end with 
our conclusions Sec. VI. 



II. THE DEBYE-HUCKEL HOLE-CAVITY 
(DHHC) THEORY REVISITED 

The one component plasma consists of N identical 
point-particles of valence v and (positive) unit charge 
q inside a volume V with a uniform neutralizing back- 
ground of charge density ~vqn B and dielectric constant 
e. As a first approximation the free energy of this system 
can be derived in the framework of the Debye-Hiickel ap- 
proach. Then, the electrostatic potential tp created by 
some ion, fixed at the origin for instance, and all its sur- 
rounding ions satisfies the spherically symmetric Poisson 
equation V 2 ip(r) = ip"(r) + *ip'(r) = -Aitp{r)/e. The 
charge density has a contribution from the central ion, 
vq5(r), a contribution from the surrounding ions which 
are distributed - within mean field theory! - according 
to the Boltzmann factor npB(r) = vqn B exp{— (3vqip(r)}, 
and finally from the charged background. Inserting this 
into the Poisson equation and linearizing the exponential 
yields the linearized Poisson-Boltzmann equation 

O A 

<tp"(r) + - 4>'{r) = n 2 i> - —vqS{r) , (1) 



where k = ^f^irtn-Q is an inverse screening length, £ — 
<?Bf 2 , £b = (3q 2 /s is the Bjerrum length, and (5 = 1/fcsT 
is the inverse thermal energy. 

The solution of Eqn. JJJ is the well known expression 
ip(r) — vqe~ Kr [er. However, the problem with Debye- 
Hiickel theory is that the condition for linearization is 
obviously not satisfied for small r, where the potential is 
large. Indeed, since all ions have the same sign of charge, 
this implies that the particle density becomes negative 
and finally diverges at the origin. This defect was over- 
come by the Debye-Hiickel- Hole theory [l^, which arti- 
ficially postulates a correlation hole of radius h around 
the central ion into which no other ions are allowed to 



penetrate. In this case the charge density is given by 

{vq (S(r) — ub) ■ r < h 
_ 4— Wj : r > h 

The hole size h is fixed by excluding particles from a 
region where their Coulomb energy is larger than k B T, 
which gives 1 + nh = (1 + 3/ii?) 1 / 3 . Once the potential at 
the position of the central ion is known, the electrostatic 
contribution to the free energy density, fpymin) , can be 
obtained by the Debye charging process [{J [l7| • 

The simple Debye-Huckel-Hole analysis of the one- 
component plasma theory offers considerable insight into 
ionic systems and is in good agreement with Monte-Carlo 
simulations when fluctuations on the charge density 
are not relevant |4^. In principle, one can attempt to 
include such fluctuations by applying the bulk density- 
functional theory in a local way. The basic idea is to 
obtain the density distribution via functional minimiza- 
tion of the free energy 

Focp[n(r)] = F PB [n(r)} + Jd 3 r f COTT (n(r)) . (3) 

The first part, the PB free energy 

F PB [n(r)] = Jd 3 r [k B Tn(r) [in (n(r)V p ) - 1] + / c i} , 

(4) 

contains the entropy of the mobile ions, the interaction 
of the small ions with the macroion potential and the 
mean-field interaction between the counterions. Here V p 
is the particle volume. The expression / corr in Eqn. © 
accounts for the correlation between the mobile ions. 
The ion distribution can be derived by minimization of 
Eqn. under the constraint of charge neutrality. Un- 
fortunately, this variational process does not lead to a 
well defined density profile if one uses /dhh(") as the 
correlation correction / corr . The reason is that /dhh(w) 
is a concave function beyond n* « 7.86/£ 3 and asymp- 
totically behaves as — n 4 / 3 . Since there is no stabilizing 
homogeneously charged background but rather a concen- 
tration of opposite charge on the macroion surface, this 
favors the development of a distribution in which all the 
ions sit on the surface of the macroion. 

The instabilities present in the DHH approach can be 
properly overcome by recognizing that the failure of this 
model is due to the too strong requirement of local charge 
neutrality imposed by the local density approximation: A 
local fluctuation leading to an increase of particle density 
implies a corresponding increase in background density. 
Therefore, the fluctuation is not suppressed by an in- 
crease in repulsive Coulomb interactions but quite on the 
contrary favored by its decrease. To circumvent the insta- 
bilities occurring at high densities, we proposed recently 
a simple solution in which one excludes the neutralizing 
background from a cavity of radius a placed around the 
central ion (for details of the derivation of the model see 
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Ref. I n this case, the charge density can be split 

in three different regions, namely 



vqS(r) 



p{r) 



-vqn B 



< r < a 
a < r < h 



en 



(5) 



V>( r ) : h < r 

4.7T 



where the hole size h is chosen such as to yield the same 
screening (i.e., the same amount of charge within h) as 
the DHH theory, which results in 



zh = [(uj - l) 3 + (na) 



31 1/3 



(6) 



with u> = (1 + 'ini) 1 ^. Using this prescription for h, the 
free energy is obtained by Debye-charging the fluid: 



DHHC 



(na)'- 



{■ 



- 2 



\ ( 7 ) 



.2(cJ 3 -l) 

ts3 



where 



fl(tD) = (E7-l)> + ^(Ej3_l) 



(8) 



Since the DHHC free energy is a convex function of 
density, /dhhc can thus be used to account for correla- 
tions within a local density approximation. 



A. The Zero Temperature Limit DHHC (0) 

The fact that the integral in Eqn. J7J) has to be solved 
numerically obstructs a direct view on how thermody- 
namic stability is actually restored. Luckily, the crucial 
point can already be seen by focusing on the limit of zero 
temperature. In this case Eqn. © gives the expression 



h = (3/47TOb 



_3\l/3 



(9) 



for the correlation hole of the DHHC theory. This conve- 
niently implies the potential to vanish outside h. In other 
words, the region a < r < h contains the right amount of 
background charge to exactly neutralize the central ion, 
and it is appropriate to refer to this limit as "complete 
screening" . The potential in the two other regions then 
simplifies considerably: 



1 + T "b-t I + ib : 0<r<a 

4>(r) = ~ A x \ 3 q ' ( 10 ) 

+ n B (l + -) - + : *<r<h 



with the dimensionless scaled density «b given by ub = 
|7ra 3 Ub- After the Debye charging process one obtains 
the following closed expression for the excess free energy 
density: 



r(0) 



riB 4a 



(11) 



Note that the limits a — > and n B — > oo do not com- 
mute: For high densities, /3/dhhc sca l es asymptotically 
like — £ne/2a, i.e., linear with density. However, in the 
limit a — * Eqn. becomes 



DHHC 



1/3 



4/3 



(12) 



and this concave scaling with density prevents it from 
being used within a local density approximation. The 
zero temperature limit thus demonstrates in a clear way 
the key role played by the cavity of size a, which excludes 
the uniform background from the vicinity of the central 



B. How to choose a proper value for a 

Before applying this strategy to various valences and 
ionic strengths, we need to specify the parameter a. If the 
counterions have a diameter d, no other charge should be 
found at a distance r < d. Therefore, in the spirit of the 
Debye-Hiickcl theory, we tentatively interpreted a in Ref. 
[l8| as the ion diameter. This choice has led to an excel- 
lent agreement with simulations when applied to rod-like 
polyelectrolytes 0] with mono-, di-, and trivalent coun- 
terions (and no added salt), but it is of course infeasible 
for point ions. In the following we suggest an alterna- 
tive way to choose a value for a which is independent of 
excluded volume arguments, and show that this choice 
yields a good description of our Monte Carlo results and 
trends. 

We already mentioned the crucial role played by a 
in maintaining the free energy convex. We also have 
seen in our discussion of the zero-temperature case that 
this is achieved because a in Eqn. © balances the 
length (3/47rriB) 1 / 3 , which is basically the mean dis- 
tance between ions. One could thus try to selfconsis- 
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tently choose a proportional to the local ion distance, 
but this would be unsuccessful: The balance would not 
work, since each density increase would shrink a propor- 
tionally, and the collapse could not be stopped. One 
thus needs a length which is local and somehow related 
to the ion density - but which does not change as the 
local ion density changes. This suggests to pick the av- 
erage distance between ions as predicted by PB theory: 
a — (3/47rnpB(^)) 1 ^ 3 - Our density functional then quite 
naturally emerges as a next order correction to the mean 
field result. 

After these general considerations on the cutoff a, let 
us continue with a practical remark. Far away from the 
charged surface the ion density is always quite low, cor- 
relations are weakly developed, and the precise value of 
a is immaterial. In fact, we only ever need a stabilizing 
cutoff close to the charged surface, where the ion density 
is largest. This suggests the following simplification: In- 
stead of using a cutoff function a(npB(f)) depending on 
the local PB density, we pick a constant a from a worst- 
case scenario, namely, the value which it has at contact. 
This then finally yields the following prescription for a: 



(47mp B (r )) 



1/3 



(13) 



In fact, since the cutoff will become important in the 
regime of strong correlations, we could even replace the 
contact density npB^o) by its limiting value 27t^bc 2 , 
where a is the density of surface charges [50|, We then 
find 



O strong coupling 



where 



8tt 2 ct 2 < 



2d 



1/3 



= 0.721 r 



-4/3 
2d 



(14) 



(15) 



is the 2d plasma coupling parameter |l5j. Formula (|14|) 
nicely demonstrates that in this limit the cavity size, 
measured in the appropriate length scale I (see also the 
scaling discussion in the Appendix), is simply another 
measure of the coupling strength. 



III. APPLICATION TO THE SPHERICAL CELL 
MODEL 

Charged spherical colloids are common and well char- 
acterizable systems for studying many electrostatic phe- 
nomena in pure culture, and they can often serve as sim- 
plified models for more complicated systems like polyelec- 
trolytes or proteins. Solutions containing such charged 
structures are indeed complicated to describe due to the 
long-range nature of the Coulombic interactions. How- 
ever, as long as these long range forces are repulsive, 
the colloids will create large correlation holes ("cells") 
around themselves which are void of other colloids. In a 



first approximation one can then decouple the macroion 
interactions and concentrate on what's going on within 
a single correlation hole an approach which is termed 
"cell model" The cell picture is known to give a 

good approximation for many realistic systems, and most 
of the physics of the system is determined by the screen- 
ing of the macroion by the microions inside a cell. As 
a test case for our theory we shall therefore consider 
a charged spherical colloid of radius ro containing Z 
charged groups, which are neutralized by point-like ions 
of valence v. This macroion is embedded in the center 
of a spherical cell of radius R, corresponding thus to a 
volume fraction <p = {ro/R) 3 of colloids. 

The thermodynamic behavior of the colloidal system 
is determined by the distribution of mobile ions around 
the macroion. This distribution is obtained by minimiza- 
tion of the free energy functional, Eqn. ©. For the 
colloidal system, the interaction of the small ions with 
the macroion and the mean-field interaction between the 
counterions are given by 



fo 



-vqn(r)(r/j(r) + ip Rx (r)^ , (16) 



where ip(r) is the total electrostatic potential at posi- 
tion r and ip& x (r) = —Zq/er is the potential due to the 
charged macroion alone. The inter-particle correlations 
are taken into account by employing / corr = /dhhc- The 
minimization itself is accomplished by numerically solv- 
ing the corresponding Euler-Lagrange-equation. Special 
care had to be taken to obtain a sufficient accuracy of 
the rapidly varying density profiles close to the colloid 
surface. 

In the following we will concentrate on two observables. 
The first is the integrated fraction of ions within a radial 
distance r from the colloid center, which is given by 



P(r) 



dr iirr 2 vqn(r) 



(17) 



We will sometimes plot P as a function of — 1/r, which is 
the Green function of the spherical Laplacian. This will 
visually expand the region close to the colloid, but it also 
has practical advantages when estimati ng t he amount of 

closely associated ions, see e.g. Ref. mo. 

Measuring all lengths in the full partition function of 
the cell model in units of I = £bv 2 reveals that the distri- 
bution function P(r) is invariant under a rescaling which 
keeps the number of counterions N = Z/v, the reduced 
colloid size r /£, and the volume fraction <fr = (r /R) 3 
constant (see Appendix) . The same holds for PB theory, 
and it is also true for DHHC theory. In the latter case 
this not only relies on the form of the DHHC free en- 
ergy correction (JJJ, but also on our particular choice of 
a. This invariance property is thus a further support for 
Eqn. 0. 

The second observable we look at is the osmotic pres- 
sure n. For PB like free energy functionals with an ad- 
ditional density term - like our / corr - it is given by |2?| 
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/3U = 



. df coir (n) 

n + n /corr(n) 

on 



i(R) 



(18) 



For the PB case, / corr = 0, this reduces to the well known 
fact that the pressure is given by the boundary density 
|53| . Since this result actually holds rigorously for the full 
restricted primitive model |50j , one could also argue that 
DHHC theory is an approximate way to calculate the 
boundary density, and then calculate the pressure from 
(3H = n(R), i.e., leave out the additional term nf — 
f. This would lead to a different result, reminding us 
that selfconsistency and consistency with other rigorous 
results cannot generally be achieved. We will always use 
the internally consistent equation l|18|) for our pressure 
calculations. 

Inserting the DHHC expression J7J) for /corr, we find 



DHHC 



(naf 



n{R) 



Alo 



257- 1 



$(57) 
fliZJ) 1 / 3 
(5J 2 - 5J)<f>(ZJ) 



(l + fi(5J)V3) (^(57)1/3 + n(57) 2 /3)2 



,(19) 



where $(57) = (57— 1) 2 + (ko) 3 lD 2 / k£. In the zero tempera- 
ture limit DHHC^ this simplifies considerably. A closed 
expression can easily be derived by combining Eqns. 
and HHJ): 



[3n°r. 



DHHC 



n(R) 



= 1 



1 + n 



-n-i/3 



1 + n" 



1 - 



1 ) 2/3 } (20) 
An 1 /3{i_3H 2 / 3 + 0(«)} ,(21) 



^ira 3 n(R). Observe that the contribution 



where n 

originating from the nf — f term is negative for all den- 
sities. 



IV. SIMULATIONAL DETAILS 

The systems we study consist of a spherical macroion 
of radius ro and (negative) central charge — Zq. Elec- 
troneutrality is ensured by the presence of N = Z/v 
point-like counterions of valence v, confined inside an im- 
permeable spherical cell of radius R. This also fixes the 
colloid volume fraction to <f> = (ro/R) 3 . No additional 
salt is added. The dielectric constant e is assumed to 
be uniform throughout the system, such that no image 
forces occur. Our choices for the system parameters 
can be found in Table [I] 

Standard canonical MC simulations following the 
Metropolis scheme ,55| were employed to sample the 
ion distributions. After an initial equilibration time of 
200 000 MC steps, where we attempted to move every 
ion once to a new position, we sampled the system for 



System 


Z 


N 


v 


r /e 


r 2 d 


npB(ro)£ 3 


a/ 1 


1 


100 


100 


1 


2 


2.5 


20.77 


0.23 


2 


120 


120 


1 


5.477 


1 


0.4090 


0.836 


3 


120 


120 


1 


2.739 


2 


8.275 


0.307 


4 


120 


120 


1 


1.826 


3 


45.07 


0.174 


5 


120 


60 


2 


1.937 


2 


7.526 


0.317 


6 


120 


40 


3 


1.581 


2 


6.976 


0.324 



TABLE I: The parameters of the simulated systems. The 
volume fraction was always chosen as <f> = (r /R) 3 = 0.8%. 




-r Q /r 

FIG. 1: Counterion distribution function P(r) for system 1 
(see Tab. GJ. The solid curve is the result of the MC sim- 
ulation, while the dash-dotted curve is the prediction from 
PB theory. The inset shows the local density n(r). The in- 
crease in the counterion condensation due to correlations is 
well captured by the DHHC theory (dashed curve) and its 
zero temperature limit DHHC' ' (dotted curve). The differ- 
ence in n(r) between the latter two is invisible on the chosen 
scale, and only DHHC is shown. 



1.3 - 2 x 10 6 MC steps, producing 1300-2000 configura- 
tions for analysis. We will measure energies in units of 
k-gT and use the coupling length i = £^v 2 as our unit 
of length (for monovalent ions under aqueous conditions 
and room temperature we would have £ = 7.14 A, and the 
unit of concentration becomes £~ 3 = 4.56 M). In the fol- 
lowing we will present MC results for the integrated ion 
distribution, Eqn. I|17ll . and for the pressure, Eqns. I|19|) 
and (J20J). 



V. COMPARISON BETWEEN SIMULATIONS 
AND DHHC THEORY 

A. Ion distribution functions 

Figure^shows the integrated charge distribution P(r) 
for system 1 from Tab. [U The solid curve is the result 
from the MC simulation, and it lies distinctly above the 
PB result (dash-dotted curve), indicating a stronger con- 



7 




1 1.25 1.5 ,2 3 4 5 

-r /r 



FIG. 2: Counterfoil distribution function P(r) for systems 
2, 3, and 4 from Tab. U] The line styles are the same as in 
Figure the counterions are monovalent, and the value of 
the plasma-parameter T2d is indicated. 



densation of ions due to correlations neglected in PB the- 
ory. Most of this enhancement of ion localization close to 
the colloid is captured by DHHC theory (dashed curve) 
or its zero temperature limit DHHC^ -'. This is also evi- 
dent from the local density n(r), which relatively to PB 
is enhanced at close proximity to the colloid, while it 
drops below PB at the outer cell boundary. From what 
we have said in section IIIII this also indicates that the 
pressure will be lower, and this is indeed what we shall 
find (see below). 

It should be noted that the 2d plasma parameter 
T2d = 2.5 is already slightly beyond the poi nt where at- 
tractions between two planes would arise [45| . We should 
not expect DHHC theory to work for significantly higher 
plasma parameters, since it cannot account for effects like 
attractions [2(|. However, we want to point out that here 
we only aim at properties of a single electrostatic double 
layer and not at phenomena arising from the interaction 
between two of them, and in fact the agreement seen in 
Fig- n is very encouraging. It is also quite pleasing that 
the significantly simpler zero temperature limit DHHC^ -* 
from Eqn. (|1 1|) yields essentially the same result as the 
full DHHC theory. 

Due to the scaling invariance of the partition function 
discussed in the Appendix, a system with e.g. divalent 
ions and Z = 200 or trivalent ions and Z — 300 (and 
properly rescaled Bjerrum lengths £b —> £b/v 2 ) shows 
exactly the same distribution function (not shown). 

In Figure [5] we show distribution functions P{r) for 
systems 2-4 from Tab.H] These have monovalent counte- 
rions and only differ in their value of the plasma param- 
eter T2d- Clearly, a larger plasma-parameter leads to an 
increased condensation (the curves are shifted up) — an 
effect which naturally is already present in PB theory. 
However, apart from this, at a larger plasma parameter 
the influence of correlations becomes more important, 




1 1.25 1.5 ,2 3 4 5 

-ro/r 



FIG. 3: Counterion distribution function P(r) for systems 
3, 5, and 6 from Tab. U The line styles are the same as in 
Figure the plasma parameter is T2d = 2, and the value of 
the counterion valence is indicated. For clarity, the PB curve 
is only shown for v = 3. 



and therefore the deviation between the PB prediction 
and the MC result increases for increasing T2d, which 
is also clearly seen in Fig. [21 Again, this effect is well 
captured by DHHC theory, which is always much closer 
to the MC data than to the PB result, even though its 
accuracy diminishes as T2d becomes large. 

In Fig. we show a "complementary" scan, in which 
we fixed the value of the plasma-parameter T2d = 2, 
but changed the counterion valence (systems 3, 5, and 
6 from Tab. P). Maybe surprisingly, an increase in va- 
lence leads to a decrease in condensation if it happens 
at constant plasma-parameter and colloid charge. If we 
had changed v from 1 to 2 and simultaneously replaced 
~~ * ^b/4 and Z — > 2Z 1 the plasma parameter would 
also have remained unchanged, but due to the scaling 
property of the partition function that would actually 
have been true for the whole distribution function. In- 
stead, we have reduced Ib — ► ^b/2 3 ^ 2 ~ ^b/2.83 (i.e., a 
little less strongly), but have failed to increase Z. The 
net result is that condensation drops slightly. However, 
since the plasma parameter, which is the best indicator 
for the strength of correlations, remains the same, the 
deviation between PB theory and MC simulation are al- 
ways about the same (not shown in the Figure). And 
as a consequence, the deviation between DHHC theory, 
which approximately accounts for correlations, and the 
MC simulation, which captures them all, is about the 
same in all three cases, and actually not very big. 



B. Osmotic Pressure 

Another strategy to check how successful our approach 
captures correlations is to compute the osmotic pressure. 
In real systems this pressure will depend on correlations 
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FIG. 4: Counterfoil density close to the cell boundary for 
system 2. The dots denote the results of the MC simulation, 
the other line styles are the same as in Fig. 



between ions of different cells, something which neither 
our theory nor actually our simulation (of a single col- 
loid!) takes into account. So in the following by "pres- 
sure" we do not, strictly speaking, refer to the bulk pres- 
sure of a colloidal suspension at some given volume frac- 
tion, but only to the pressure exerted on the rigid wall 
at r = R of our cell model. 

Within the simulations, the pressure is given by the 
contact density at r = R, which was obtained by fit- 
ting the MC density profile close to the cell boundary 
to a quadratic expression n(r) = c\ + ci{r — R) 2 . An 
example for how the simulated densities compare to the 
PB approximation, our analytic DHHC approach, and its 
simpler zero temperature limit DHHC^ -*, can be found 
in Fig. H 

Table [H] shows the predictions for the pressure in 
the case of point-like ions given by PB-, DHHC-, and 
DHHC -theory, as well as by MC simulations. In the 
case of PB-theory and MC the pressure is simply the 
density at the outer boundary, while for the DHHC ap- 
proach we employ Eqn. (JTSJ and for DHHC (0) Eqn. i|2Tjl. 
These data, as well as Fig. demonstrate that - as an- 
ticipated - the simulated pressures lie below the PB pre- 
diction. This decrease in pressure is rather accurately 
captured by our functional /dhhc and by its zero tem- 
perature limit, Jd°h HC , with the MC result lying signifi- 
cantly below PB and (in these cases) below DHHC and 
above DHHC^. The difference between the two corre- 
lation corrected approaches is consistent with the idea 

that entropic effects neglected in /q°hhc wou ld push ions 
away from the macroions, or in other words, that the 
zero temperature limit implies stronger correlations than 
the regular DHHC theory and therefore yields even lower 
pressures. 



VI. CONCLUSIONS 

In this paper we showed how to apply our previously 
proposed local density functional approach based on a 
stable correlation correction to a spherical macroion con- 
fined in a spherical cell. One of the crucial parameters 
in this theory is the size a of the exclusion cavity of the 
background charge density. For point-like ions, we sug- 
gest to associate the exclusion region with the mean dis- 
tance between ions as predicted by PB theory, and for 
simplicity use the value present at colloidal contact. 

By going to the zero temperature limit we were able 
to derive an even simpler free energy functional -Fdhhc > 
which is almost as good as the full DHHC theory, but 
much easier to handle. We also derived exact expres- 
sions for the osmotic pressure in this system. We suc- 
cessfully compared our predictions to simulations of the 
same model and compared the integrated counterion den- 
sity and the osmotic pressure values for two complemen- 
tary "scans" of the coupling strength, namely valence 
and plasma parameter. We demonstrated that our local 
density functional approach based on a stable correla- 
tion correction leads to a major improvement over the 
PB prediction. 
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APPENDIX 

The canonical partition function Z of the colloid sur- 
rounded by its counterion is given by: 

i— 1 

where N = Z/v is the total number of counterions and 
the Hamiltonian TL = T + V splits into kinetic and poten- 
tial degrees of freedom. In the classical description em- 
ployed here the kinetic part T will contribute the usual 
factor A~ 3N to the partition function, where A is the ther- 
mal de Broglie wavelength. The potential energy can be 
expressed as 

v --"?R + 5gi^T- < 23 > 
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Sys. 


r 

/SIIpb^ 




3 

/OTdhhc^ 


3 

/3n DHH c^ 


„— 

(3Hmc£ 




ITpb/IImc 


IlDHHc/nMC 


ttO /tt 
llDHHC/iiMC 


1 


2.98 x 10" 


-3 


2.56 x 10~ 3 


2.39 x 10~ 3 


2.53(3) x 10" 


- 3 


1.18 


1.01 


0.94 


2 


3.74 x 10" 


4 


3.58 x 10~ 4 


3.44 x 10~ 4 


3.55(4) x 10" 


-4 


1.05 


1.01 


0.97 


3 


1.58 x 10" 


•3 


1.42 x 10~ 3 


1.33 x 10~ 3 


1.38(3) x 10" 


-3 


1.14 


1.03 


0.96 


4 


3.61 x 10" 


-3 


3.02 x 10" 3 


2.81 x 10~ 3 


2.95(5) x 10" 


-3 


1.22 


1.02 


0.95 


5 


3.09 x 10" 


•3 


2.72 x 10" 3 


2.53 x 10" 3 


2.63(6) x 10" 


-3 


1.17 


1.03 


0.96 


6 


4.10 x 10" 


•3 


3.96 x 10" 3 


3.71 x 10~ 3 


3.83(9) x 10" 


-3 


1.07 


1.03 


0.97 



TABLE II: The values of the various pressures (in units feT /£ 3 ) for the systems 1-6. The MC errors have been conservatively 
estimated from the fluctuations of the measured density around the fit close to the cell boundary. The last three columns 
display the ratio between the theoretical and MC values, illustrating which theories over- or underestimate the pressure, and 
by how much. 



After rescaling all length by £, i.e. introducing x := r/l, 
the total partition function can be written as 



1 / f\ 3N /• a; o/0 1/3 

Jx « k 



cxpi-A^y -L + iy 1 



(24) 



In this form it becomes evident that appropriately scaled 
thermal observables like the integrated charge density 



(measured in units of £~ 3 ) or the pressure (measured 
in units of k^T£~ 3 ) are invariant under system changes 
which leave the number of counterions N, the rescaled 
colloid size x = r /£, and the volume fraction <f) fixed. 



Poisson-Boltzmann theory shows the same invariance 
property, as does the approximate density functional the- 
ory we are proposing in this paper. 



[1] G. L. Gouy, J. de Phys. 9, 457 (1910). 

[2] D. L. Chapman, London and Edinburg Phil. Mag. and J. 

of Sci. 25, 475 (1913). 
[3] B. V. Derjaguin and L. D. Landau, Acta Physicochim. 

(USSR) 14, 633 (1941). 
[4] E. J. Verwey and J. T. G. Overbeek, Theory of the stabil- 
ity of Lyophobic Colloids (Elsevier, Amsterdam, 1948). 
[5] P. Debye and E. Hiickel, Phys. Z. 24, 185 (1923). 
[6] L. Belloni, Phys. Rev. Lett. 57, 2026 (1986). 
[7] J. P. Hansen and I. McDonald, Theory of Simple Liquids 

(Academic, London, 1990). 
[8] R. D. Groot, J. Chem. Phys. 94, 5083 (1991). 
[9] M. Lozada-Cassou and J. E. S. Sanchez, Chem. Phys. 
Lett. 190, 202 (1992). 
[10] R. R. Netz and H. Orland, Eur. Phys. J. E 1, 67 (2001). 
[11] H. Lowen, J. -P. Hansen, and P. A. Madden, J. Chem. 

Phys. 98, 3275 (1993). 
[12] S. Nordholm, Chem. Phys. Lett. 105, 302 (1984). 
[13] E. E. Salpeter, Ann. Phys. 5, 183 (1958). 
[14] R. Abe, Prog. Theory Phys. 22, 213 (1959). 
[15] M. Baus and J.-P. Hansen, Phys. Rep. 59, 1 (1980). 
[16] R. D. Groot, J. Chem. Phys. 95, 9191 (1991). 
[17] R. Penfold, S. Nordholm, B. Jonsson, and C. E. Wood- 
ward, J. Chem. Phys. 92, 1915 (1990). 
[18] M. C. Barbosa, M. Deserno, and C. Holm, Europhys. 

Lett. 52, 80 (2000). 
[19] N. Gr0nbech- Jensen, K. M. Beardmore, and P. Pincus, 

Physica 261A, 74 (1998). 
[20] E. Allahyarov, I. D'Amico, and H. Lowen, Phys. Rev. 

Lett. 81, 1334 (1998). 
[21] Y. Levin, J. J. Arenzon, and J. F. Stilck, Phys. Rev. Lett. 
83, 2680 (1999). 



[22] R. Messina, C. Holm, and K. Krcmer, Phys. Rev. Lett. 
85, 872 (2000). 

[23] R. Messina, C. Holm, and K. Kremer, Europhys. Lett. 

51, 461 (2000). 
[24] J. Neu, Phys. Rev. Lett. 82, 1072 (1999). 
[25] J. Sader and D. Y. Chan, J. Colloid Interface Sci. 213, 

268 (1999). 

[26] E. Trizac, Phys. Rev. E 62, R1465 (2000). 
[27] G. Tellez and E. Trizac, J. Chem. Phys. 118, 3362 (2003). 
[28] R. van Roij and J.-P. Hansen, Phys. Rev. Lett. 79, 3082 
(1997). 

[29] P. B. Warren, J. Chem. Phys. 112, 4683 (2000). 

[30] J. F. Dufreche, T. O. White, and J.-P. Hansen, Mol. 
Phys. 101, 1741 (2003). 

[31] H.-H. von Griinberg, R. van Roij, and G. Klein, Euro- 
phys. Lett. 55, 580 (2001). 

[32] A. Diehl, M. C. Barbosa, and Y. Levin, Europhys. Lett. 
53, 86 (2001). 

[33] M. Deserno and H.-H. von Griinberg, Phys. Rev. E 66, 
011401 (2002). 

[34] M. N. Tamashiro and H. Schiessel, J. Chem. Phys. 119, 
1855 (2003). 

[35] F. Oosawa, Poly electrolytes (Marcel Dekker, New York, 
1971). 

[36] G. Manning, J. Chem. Phys. 51, 924 (1969). 
[37] M. Le Bret and B. Zimm, Biopolymers 23, 271 (1984). 
[38] M. Deserno, C. Holm, and S. May, Macromoleculcs 33, 
199 (2000). 

[39] S. Alexander et al, J. Chem. Phys. 80, 5776 (1984). 
[40] L. Belloni, Colloids and Surfaces A 140, 227 (1998). 
[41] E. Trizac, L. Bocquet, and M. Aubouy, Phys. Lett. Lett. 
89, 248301 (2002). 



10 



[42] E. Trizac, M. Aubouy, and L. Bocquet, J. Phys. Cond. 

Mat. 15, S291 (2003). 
[43] E. Trizac, L. Bocquet, M. Aubouy, and H.-H. von 

Griinberg, Langmuir 19, 4027 (2003). 
[44] M. Aubouy, E. Trizac, and L. Bocquet, J. Phys. A: Math. 

Gen. 36, 5835 (2003). 
[45] A. G. Moreira and R. R. Netz, Eur. Phys. J. E 8, 33 

(2002). 

[46] N. Gr0nbech- Jensen, R. J. Mashl, R. F. Bruinsma, and 
W. M. Gelbart, Phys. Rev. Lett. 78, 2477 (1997). 

[47] M. Deserno, A. Arnold, and C. Holm, Macromolecules 
36, 249 (2003). 

[48] S. G. Brush, H. L. Sahlin, and E. Teller, J. Chem. Phys. 



45, 2102 (1966). 
[49] M. N. Tamashiro, Y. Levin, and M. C. Barbosa, Physica 

A 268, 24 (1999). 
[50] H. Wennerstrom, B. Jonsson, and P. Linse, J. Chem. 

Phys. 76, 4665 (1982). 
[51] A. Katchalsky, Pure Appl. Chem. 26, 327 (1971). 
[52] L. Belloni, M. Drifford, and P. Turq, Chem. Phys. 83, 

147 (1984). 

[53] R. A. Marcus, J. Chem. Phys. 23, 1057 (1955). 
[54] R. Messina, J. Chem. Phys. 117, 11062 (2002). 
[55] N. Metropolis et ai, J. Chem. Phys. 21, 1087 (1953). 



